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ABSTRACT 

In 1961, Sperling linearized and regularized the differential equations of motion 
of the two-body problem by changing the independent variable from time to fictitious 

time by Sundman’s transformation (r = —) and by embedding the two-body energy 

integral and the Laplace vector. In 1968, Burdet developed a perturbation theory 
which was uniformly valid for all types of orbits using a variation of parameters 
approach on the elements which appeared in Sperling’s equations for die two-body 
solution. In 1973, Bond and Hanssen improved Burdet’s set of differential equations 
by embedding the total energy (which is a constant when the potential function is 
explicitly dependent upon time.) The Jacobian constant was used as an element to 
replace the total energy in a reformulation of the differential equations of motion. In 
the process, another element which is proportional to a component of the angular 
momentum was introduced. 


Recently trajectories computed during numerical studies of atmospheric entry from cir- 
cular orbits and low thrust beginning in near-circular orbits exhibited numerical insta- 
bility when solved by the method of Bond and Gottlieb (1989) for long time intervals. 
It was found that this instability was due to secular terms which appear on the right- 
hand sides of the differential equations of some of the elements. In this paper, this 
instability is removed by the introduction of another vector integral called the delta 
integral (which replaces the Laplace Vector) and another scalar integral which remove 
the secular terms. The introduction of these integrals requires a new derivation of the 
differential equations for most of the elements. For this rederivation, the Lagrange 
method of variation of parameters is used making the development more concise. 
Numerical examples of this improvement will be presented. 


This work was performed for NASA-JSC Houston, Texas under Contract No. NAS9- 
17885. 
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1.0 Summary 


In 1961 Sperling linearized and regularized the differential equations of motion of the two-body prob- 
lem by changing the independent variable from time to fictitious time by Sundman’s transformation 

(r = — ) and by embedding the two-body energy integral and the Laplace vector which is also an 

integral of the motion into the Newtonian form of the differential equations of motion. The solution of 
Sperling’s differential equations was uniformly valid for all types of orbits. In 1968, Burdet developed 
a perturbation theory using a variation of parameters approach on the 14 elements which appeared in 
the two-body solution. In 1973, Bond and Hanssen improved Burdet’s set of differential equations by 
using the total energy of the perturbed system as a parameter instead of the two-body energy and by 
reducing the number of elements to 13. In 1989 Bond and Gottlieb embedded the Jacobian integral, 
which is a constant when the potential function is explicitly dependent upon time as well as position in 
the Newtonian equations. The Jacobian constant was used as an element to replace the total energy in 
a reformulation of the differential equations of motion. In this process, another element which is pro- 
portional to a component of the angular momentum is introduced. This brought the total number of 
elements back to 14. In this paper the Laplace vector is replaced by another vector integral as well as 
another scalar integral which remove small secular terms which appear in the differential equations for 
some of the elements. 

2.0 Introduction 


The non-linear differential equations of motion for the cartesian coordinates of the two-body problem 
can be regularized and linearized by the three-step procedure of changing the independent variable form 
time (t) to fictitious time (s) by the application of the Sundman transformation, embedding the Laplace 
integral and embedding the Jacobian integral. 


By regularization we mean the removal of all singularities, and by linearization we mean that the 
differential equations for the cartesian coordinates are transformed to harmonic oscillators. Previously, 
regularization and linearization were done by Burdet (1968) by embedding the two-body energy which 
is constant only for the two-body problem and by Bond and Hanssen (1973) by embedding the total 
energy which is a constant when the two-body system is perturbed by a conservative potential (function 
of position only). In Bond and Gottlieb (1989), the Jacobian integral, which is a constant for the case 
of the two-body system perturbed by a potential function that is explicitly dependent on time as well as 
position, was embedded in the Newtonian equations. All three of these approaches reduce to the same 
system of equations in the absence of perturbations. 


Recent numerical studies on atmospheric entry from near circular orbits and on low thrust in near circu- 
lar orbits exhibit numerical instability when solved by the method of Bond and Gottlieb (1989) for long 
time intervals. These two cases are similar since both have persistent, tangential, non-conservative per- 
turbations. It was found that this instability was due to secular terms which appear on the right hand 
sides of the differential equations of some of the elements. In this paper this instability is removed by 
the introduction of another vector integral of the motion and another scalar integral which remove the 
secular terms. The introduction of these integrals which were included by Burdet (1968) require a new 
derivation of the differential equations for most of the elements. For this rederivation the Lagrange 
method of variation of parameters is used making the development more concise. 

2.1 The Differential Equations of Motion in the Fictitious Time 


The differential equation for perturbed two-body motion is 


r + 


5i-2 


( 2 . 1 ) 


where r_ is the position vector of one of the masses with respect to the other in cartesian coordinates 

and r is the magnitude of r and ( ) = Also the gravitational constant is 

at 
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\a = G{M+m) ( 2 - 2 ) 

where G is the universal gravitational constant and M and m are the masses of the two bodies. The 
quantity F is the perturbation which can be expressed by, 

F=P-^~V(r,t) < 2 - 3 ) 

— ” dr 

where V( r, t) is the potential due to perturbing masses and P is any perturbative acceleration which is 
not derived from a potential. 

Equation (2.1) can be linearized (except for the perturbation) in three steps: 

STEP (1) Change the independent variable from time (t) to fictitious time (s) according to the transfor- 
mation 



ii 

43|43 

(2.4) 

The derivatives of r with respect to t become 

r = r / r 

(2.5) 

where ( ) - — p-, and 
as 



r = r" / r 2 - rV / r 3 

(2.6) 

where 

r = r ■ t / r 

(2.7) 


STEP (2) Embed the integral called the Laplace vector (a constant when F= 0) 




rtr 


which becomes 




when the new independent variable s is used. 

STEP (3) Embed the energy integral (a constant when F=0) 

_ 2 y. 


a* = 


r • r 


which becomes 


2|i. 1 - - 


( 2 . 8 ) 


(2.9) 


( 2 . 10 ) 


( 2 . 11 ) 


when the new independent variable is used. Note that 

a k =-2h k ( 2 - 12 ) 

where h k is the two-body or Keplerian energy. Using these three steps in order, equation (2.1) becomes 
r" + a k r = - \i£ + r 2 F ( 2 - 13 ) 

which is the differential equation for the position vector r. By taking the dot product of equation 
(2.13) with the position vector r we obtain 

t" + a k r = p + r r • F (2.14) 

which is the differential equation for the distance r. We now change from the energy integral a k to 
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(2.15) 


the Jacobian integral a j (Bond and Gottlieb (1989)) which is given by 

a.j = ct k + 2a - 2 V(r, t) 

where a is called the axial element and is defined by 


CT = 0 >-(rxr) (2.16) 

The vector co is the constant rotational rate of the central attracting body or orbital rate of a third body 
giving rise to the perturbing potential V( r, t). In Section 4.0 it will be shown that a, » constant 
when /£/ = 0 and that a = constant when /oV = 0. Solving equation (2.15) for a* and substituting into 
equations (2.13) and (2.14) we obtain 


and 


r + a.] i = - |i£ + r 2 F + 2(a - V( r, t))r • - p£ + Q 


(2.17) 


r" + a,r = \i + rr ■ F + 2(a - V( r , r))r * p + jQ • r (2.18) 

Note that all of the perturbation terms have been moved to the right side in equation (2.17) and (2.18). 


Equation (2.17) and (2.18) are coupled only through the perturbation terms. We will refer to equation 
(2.17) as the spatial differential equation since its solution provides position and velocity. We will refer 
to equation (2.18) along with equation (2.4) as the temporal differential equations since their solutions 
provide time. Note that when there are no perturbations (that is IF I = 0 and /oV = 0) then ve have the 
two-body differential equations 


and 


L + a tL - ~ 


(2.19) 


r" + a jr = p 

and the Jacobi constant and Keplerian energy become the same 


( 2 . 20 ) 


«/ - a* 


3.0 Two Body Solution 

The differential equation of motion for the two-body problem in the fictitious time was shown in the 
previous section to be 

r' + a y r=-p£ (3.1) 

The solution of (3.1) in terms of the Stumpff functions of Appendix B is 

r=r p c 0 + r^sc , - \lts 2 c 2 ( 3 . 2 ) 

where and are the initial values of r and r_ , and the Stumpff functions are c t = c t (ajs 2 ). This 
can be verified by direct substitution of (3.2) into (3.1) and using the derivatives of the Stumpff func- 
tions 

C 0 =- OtySC] 

sc[+ci=c 0 ( 3 . 3 ) 

SC 2 + 2c 2 = Cl 

The first derivative of (3.2) which is the "velocity" in the fictitious time is 

1 = - (a/£«, + P§)ic, + £p'c 0 (3.4) 

In place of pe which is a constant of the motion we define the constant "delta vector" 
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(3.5) 


5 = - ctyr, - |i£ 

Now using the Stumpff function identity 
C g + OLjS 2 C 2 ■ 1 

and equation (3.5) and (3.2) we obtain 

l=L> + b,sc x + 8j 2 c 2 


(3.6) 

(3.7) 


similarly equation (3.4) becomes 

r = £,'c 0 + 8iCi 

The differential equation of motion for the distance r was shown in the previous section to be 


r" + a ]T = |A 

The solution of equation (3.9) is similar to that for (3.1). In terms of Stumpff 

T = TgCg + r'gSC J + ^ ^ 2 


(3.9) 

functions the distance is 
( 3 . 10 ) 


and its derivative is 


T - (M- ~ r o aj)SCi + fgCg 


(3.11) 


Now define the constant 

y = \i-r 0 aj 

which we substitute for (j. in equation (3.10) along with the identity of equation (3.6) to obtain 

r =r 0 + r' 0 sc i + ys 2 c 2 


(3.12) 

(3.13) 


Similarly equation (3.1) becomes 

r = T’ 0 Cg + ysc\ ( 3 - 14 ) 

Now substitute equation (3.13) for r in the independent variable transformation, equation (2.4), to 
obtain 


dt =r 0 ds + sc ids + ys 2 c 2 ds 
Now use the integration formula 

J s m c m ds = s' n+1 c m+1 

which is from Appendix B to obtain the equation for time (Kepler’s equation), 


t =t 0 +r 0 s +r 0 s*c 2 + ys c 3 


(3.15) 


(3.16) 


where t 0 is the initial time. 

The integration constants which were introduced in this section are r^, r^,, r g , r 0 , t 0 . The new constant 
8 simply replaces the Laplace vector which is a constant of two-body motion through the definition 
( 3 . 5 ). Similarly we note that the constant y replaces the gravitational constant (equation (3.12)). The 
introduction of the constants 8 and y was done by Burdet (1968). This fact was noted by Bond and 
Hanssen (1973). The Jacobian element ot j is the same as the two-body energy parameter a* in the 
unperturbed case is also a constant of the motion. In addition we have the axial element cr which is 
also a constant of the motion (see equation 2.16). This is a total of 15 constants of the motion. 


The constants and 8 will be treated as orbital elements associated with the spatial differential 

equation (2.17) and r 0 ,F o ,y,t 0 will be treated as orbital elements associated with the temporal 
differential equations (2.4) and (2.18). 
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4,0 The Differential Equations For The Elements 


When perturbations are present the elements are no longer constant. First we derive the differential 
equation for the axial element a. Differentiate equation (2.16) with respect to time and substitute equa- 
tion ( 2 . 1 ) and (2.3) to obtain 


a = o)-(r xf') = a- rxf = rr x 
now use equation (2.4) to change to fictitious time 


P - ^V(r, t) 


a = r© • r x 


J-Jrva.o 


(4.1) 


(4.2) 


Clearly <y — constant when /to/ — 0. Now we derive the differential equation for the Jacobian element 
d j. Differentiate equation (2.15) with respect to time to obtain 


dy = a* + 2 df - 2 
From equations (2.10) and (2.1) 




dr 


a* = - 2r 


P~f v (r,t) 


and from Bond and Mulcihy (1988) also Bond and Gottlieb (1989) 

^■V(r,t)=-(orx^-V(r,t) 

and from equation (4.1) the expression for dy becomes 


dy = 2(-r + a x r) ■ P (4.3) 

Now use equation (2.4) to change to fictitious time 

ay = 2(-r' + r w x r) • P ( 4 . 4 ) 

Note that ay = constant when IPI = 0. The Jacobian constant ay will be treated as an orbital element 
for both the spatial and temporal equations since a y appears in the two-body equations (2.19) and 
(2.20). Even though we have already obtained the differential equation for a y (equation (4.4)) we must 
include it in the variation of parameters procedures of the spatial and temporal equations. The axial 
element <y appears only as a perturbation in equations (2.17) and (2.18). We have also obtained the 
differential equation for a (equation (4.2)). We will include a in the variation of parameters procedure 
for convenience and completeness. 


Even though the Laplace vector will be eliminated as an element we will need the derivative of the 
Laplace vector as a perturbation. This derivative as found by differentiating equation (2.8) will respect 
to time, then using equation ( 2 . 1 ) to eliminate r, and finally using equation ( 2 . 4 ) to obtain 

= 2 (r' • F)r - (r • F)r' - (r ■ r')F (4.4a) 


4.1 Spatial Elements 

Now we apply the variation of parameters method of Lagrange to equations (2.17), (4.2) and ( 4 . 4 ). 
Define 


~L 

* 2=1 

£3 = - ay r - H£ = -Xjjc, - p£ (4.5) 
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X 4 = O 
*5 = a J 

Now differentiate equations (4.5) and use (2.17), (4.2) and (4.4) to obtain 

£l = £2 

£2 = £3 + 2 = £3 + Q.z 
£3 = ~ x 5£2 - («/£ + M£) = - *5£2 + <£3 
X 4 = o' = G 4 
*5 = <*/ = G 5 

Where Gi = 0. Equations (4.6) can be separated into unperturbed (i.e., two-body or Keplerian) and per- 
turbed parts, that is into the form of x ' - F + G , making them suitable for Lagrange’s variation of 
narameters method as given by Appendix A. In this form equation (4.6) becomes 


(4.6) 


F 1 
x[ 


p ^ 

£2 


r 

*2 


£3 


g 2 

£3 

= 

~* 5£2 

+ 

G 3 

X 4 


0 


G4 

X5 
L. J 


0 

L J 


G 5 


(4.7) 


where G lt G 2 , G 3 , G 4 , G 5 are defined from equations (4.6). The array of constants, which will become 
the new dependent variables, is defined as 

c T = (a T , g r , 8 T , a, a y ) ( 4 - 8 ) 

where 

a = £p = £i(0) 

e=^'=£ 2 ( 0 ) ( 4 - 9 > 

8 = - ay a - |i£ = x 3 ( 0 ) 

and of course a and a y which have already been established as constants of the motion. The 
differential equation for c , has the form, — c = G where 


dx 

dc 


3*i 

3xi 

3xi 

axr 

3X! 

3a 

dg 

d8 

da 

day 

3x 2 

3x 2 

3x 2 

dx 2 

3x 2 

da 

dl 

d8 

3a 

day 

3x3 

dx 3 

dx 3 

dx 3 

dx 3 

da 

d§ 

ds 

3 a 

day 

dx 4 

3x 4 

dx 4 

3 X 4 

dx 4 

da 

« 

d8 

3a 

3a 2 

dx 5 

3x 5 

dx 5 

3x5 

3x 5 

3a 


d8 

3a 

day 


( 4 - 10 ) 


Noting from Section (3.0) that 


£1 = r = a + gsci + 8r 2 c 2 
£ 2 = r = £c 0 + 8rci r 


( 4 -ll) 
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and from equation (3.5), (4.5) and (4.9) 

x 3 = a y (g-r) + 5 

also 

x 4 = a 

*5 = 0/ 

The differential equations become 

I Isci Is 2 c 2 0 

[0] ic 0 i sc , o 

[0] -Ia.jsc x Ic 0 0 


(4.12) 


dr 


da.j 

3/ 


da.j 

dx 3 


o T o T 

o T o r 


da.j 

o T i o 

(f o 1 


a 


G 1' 

e' 


G 2 

8' 

= 

G 3 

a 


g 4 

«/ 


g 5 


(4.13) 


where we have used the identity from Appendix B 

c 0 = 1 - a js 2 c 2 

also, I is the 3 by 3 identity matrix; [0] is the 3 by 3 null matrix; 0 is a column vector with 3 com- 
ponents; 0 is a row vector with 3 components. Equation (4.13) yields the equations 

a + gsci + 8 j 2 c2 + = 0 

OCCy 


j3'c 0 + 5'sci + <ij~~ = £ 

dCty 

- JJ'a ijsci + 5'c<, + a'j ^“ 3 


(4.14) 


= - a } r_ - p£ 


9a/ 

a = rm - rx F 

ctj = 2(-/ + rox r ) ■ P 

where we have restored the original notations for G x , G 2 , G 2 , G 4 , G5. Now using the partial deriva- 
tives. 


dr _ 9 c 1 

9a/ ^ 9a/ 


+ 


+ 8 s 


- 9a y 
9c 1 


-9a/ 

dr 


9r' _ „ dc c 

9a/ ^ 9a y 

dx 2 

-5 — = a - r - a, - — 

9a/ — - day 

where the Stumpff function derivatives are 

1 2 

T — = - — J Cl 

9a/ 2 


(4.15) 


(4.15a) 


1 


9a/ 2a/ 


(c*_! - kc k ) , £>1 


and other Stumpff function identities from Appendix B equation (4.14) can be solved simultaneously, 
omitting several algebraic steps to give 
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(4.16) 


a = - Qsci - ii£s 2 c 2 - aj |ar 2 c 2 + 2]}.s 3 e 3 + y8j 4 c 2 

§ = Qpo + M£ ^ci + a'j |ajc , + - Ss 3 (22 3 - c^J 

8 = Qajsc x - |j £c 0 + ajf- ac 0 + 2 o,£j 3 ? 3 + y8a/j 4 c 2 


a = r a> • r x F 

Oj = 2( - r ' + rco x r ) • P 

where £j = cj(4ays 2 ) as discussed in Appendix B. In the reference Bond and Gottlieb (1989) the 
coefficient of the factor aja in the differential equation for § had a secular term. This term does not 
appear in equation (4.16). Note that the Laplace vector (p t) has been entirely removed from the formu- 
lation. The derivative of the Laplace vector (pc ) remains but only as an abbreviation for the perturba- 
tions given in equation (4.4a). 


4.2 Temporal Elements 

Now we apply Lagrange’s variation of parameters method to equations (2.18), (2.4) and (4.4). Define 

y i = r 

yi=r 

y 3 = H-(Xyr (4.17) 

y 4 = « 
y s = a, 

Note that a.j is the only element which is common to both the spatial and temporal systems. Now 
differentiate equations (4.17) equation (2.18), (2.4), and (4.4) become 

y'i =yi 

yt =73+ -Q. L = )'3 + «2 

y'i =-ysy2-ajr = -y 5 y2 + g3 ( 418 ) 

y 4 =>-1 

ys =aj = g s 

Where = 0 and g 4 = 0. Equation (4.18) can also be expressed in the form y = f + g 


y{ 


’ i 

y 2 


r % 

g 1 

y'i 




8 2 


= 

~y $y 2 

+ 

g 3 



y i 


8 4 

y s 


0 

L J 


8 5 


where g t , g 2 , g 3, g 4 , g 5 are defined by equation (4.18). The array of constants which will become the 
new dependent variables are 

K r = (fl,fi,Y, x, aj) (4.20) 

where 


a = r 0 =yi( 0 ) 
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(4.21) 


b =r' = y 2 ( 0 ) 

Y=M - a ja = y 3 (0) 
x = t 0 = y 4 (.0) 

and ay has already been established as a constant of the motion. 

The differential equations for k (having the form -^--k = g ) becomes 

0K 


dyi 

3y, 

dy i 

3yi 

3yj 

da 

db 

dy 

3r 

9a y 

dy 2 

dy 2 

dy 2 

dy 2 

dy 2 

da 

db 

dy 

3r 

3ay 


dy 3 

dy 3 

dyj 


da 

db 

dy 

9t 

9ay 

dy 4 

dy 4 

dy 4 

d>4 

3y 4 

da 

db 

dy 

3t 

day 

dys 

dy 5 


3y 5 

3y 5 


| da db dy dx da j 

but from equations (3.12), (3.13), (3.14), (3.15) and (4.21) 
yi = r = a + bsc t + ys 2 c 2 
y 2 = r' = bc 0 +ysc { 

y 3 = IX- cijr = y + aja - a 2 r = y+ctj(a - r) 

>4 = t = x + as + bs 2 c 2 + ys 3 c 3 

y 5 = ay 


So we can evaluate the matrix elements in (4.22) to obtain 

1 sci s 2 c 2 0 


O n * 

c 0 sc x 0 - — a 8 1 

oOCy , 

b 82 

dy 3 

0 -ajsci c 0 0 Y = £3 

2 3 . dt Z . gA 

5 sc * s C3 1 ^7 H N 


Lo o o o i J 

Equation (4.24) when expanded yields, 

Q + i 5Ci + Y^C? + Ot; "r = 0 

oa y 

b c 0 + y sc i + aJ^—~ = —Q • t 
day r ~ 

* 'ty 3 

- b a } sc x + y c 0 + ay-r — = - ray 
day 


(4.22) 


(4.23) 


(4.24) 


(4.25) 
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as + b s 2 c 2 + ys c 3 + x + a J~^~ = 0 

ay = 2( - / + r (0 x r) ■ P 

Where we have restored the original notations for g 2 , 83 8 s- We evaluate the partial derivatives in 

equations (4.25) using equations (4.23) 

dr , &i ,dc 2 

t , — = bs ^: — + ys 2 - — 
day day day 


dr' , 3c 0 3c i 

— = b - — + ys-r — 
day 3ay 3ay 

dy 3 dr 

= a - r - a.] 


(4.26) 


day 

dr 


day 


3c 


dc - 


, OC 2 3 0< - 3 

, = in j — + Y S t — 

day day day 

where the Stumpff function derivatives are given by equations (4.15a). Equations (4.25) can be solved 
simultaneously for the derivatives, 


a = - — r ■ Qsc i - ay 


1 


as l c 2 + 2 bs 3 l-i + -rys 4 c 2 


b = —r ■ Qc 0 + a j asc : + bs 2 £ 2 - ys 3 (2c 3-C1C2) 
r~ L -i 


(4.27) 


y = — r • gay sc, + a y 


- ac 0 + 2bajs 3 d 3 + ^-yajs 4 ci 


x = —r ■ Qs 2 c 2 + a y 


as 3 c 3 + ^-fe 4 C2 - 2ys 5 (c 5 -4c5) 


As in the development of equations (4.16) the Stumpff function identities of Appendix B have been 
used. In the reference Bond and Gottlieb (1989) the coefficient of the factor ay a in the differential 
equation for b had a secular term. This term does not appear in equation (4.27). 

It is useful to note that 

H = y + a ja 

is an integral of the system of equations (4.27). From equations (4.27) it is easy to show that 
y + aa.y + art] = 0 
which can be integrated to give 

y + a ay = constant (4.30) 

By comparison of equation (4.30) to equation (4.21) the constant of integration is the gravitational con- 
stant p. Therefore it is not necessary to compute y from its differential equation. We can compute y 
from equation (4.28), 

Y=p-a y a (4.31) 

5.0 Minimization Of Perturbations 


(4.28) 

(4.29) 


The variation of parameters approach is not dependent on the magnitude of the perturbation. No 
assumption on the size of the perturbation is required in order that the method be rigorous. However, 
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small perturbations enhance the efficiency, speed, and accuracy of any perturbation method. In the 
method described in this paper, the embedding of the Jacobi integral has the effect of introducing a per- 
turbation parameter /co/ that is the rotational speed of the planet, or the mean motion of the perturbing 
third body. To prevent this perturbation from becoming too large the following modification is offered: 

Let, 


c = c 0 + Act (5.1) 

where a 0 is the initial value of a and Act is the change in a. In effect we can let Ac replace a so that 
the differential equations reflect only changes in a. Substitute equation (5.1) into equation (2.17) to 
obtain 

r" + ajr = - p£ + r 2 F + 2 (a 0 + Ac - V(r , t))r 
Now since a* is constant we can move it to the left side of this differential equation to get 

r" + (a, - 2a* )r = - |i£ + r 2 F + 2(Aa - ^(r, t))r (5.2) 

Similarly equation (2.18) becomes 

r" + (a j - 2 a Q )r = p + rr ■ F + 2(Aa - V(r, ())r (5.3) 

This change does not affect the outcome of the variation of parameters approach taken here. This 

change is only a computational convenience and is in effect in the computational procedure of Section 
6.1 where the element ctj is actually aj - 2a 0 and a is actually Aa. Note that the initial value of Aa 
is 

Aa = 0 (5.4) 


6.0 Application 

In this section the most important equations are collected and listed in a logical order suitable for com- 
putation, Also two numerical examples are presented, 

6.1 Computational Procedure 
Given r,, t 0 find r(t) and v(t). 

STEP 1 Initialization 


5=0 

r o =(L> ' L> ) 1/2 
a = r 0 

*> = L> ■ v* 

£ = « Y* 


Evaluate Perturbations V 0 , 


dV 

*J. 





-2 V, 
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5 = - (v 0 • Yp)!j> + <L> ■ + rL> ~ «/£p 

'O 

CT = 0 

STEP 2 Transform Elements to Coordinates 

r = a + gjc i + Ss 2 Ci 
r = §c„ + Sic i 
£3 = 0 / (a - r) + 5 
7= 4 -a y a 
r = a + foci + yr 2 C2 
v = r'/r 
r = bc 9 + ysci 
t = T + as + 6j 2 C2 + 7^ 3 C3 

STEP 3 Evaluate Differential Equations For The Elements 


F = 


£-f 


2 = + 2r(-V + a ) 

aj 

Mi' 


2( - r ' + r <0 x r ) • P 

2(r' F)r-(r- F)/ - (r • r')F 


a = - Qsc 1 - M£j 2 c 2 ~ 0 j 


as 2 c 2 + 2§s 3 <J 3 + 4-§r 4 C2 


j}' = QPo + M£«i + aj |jxjci + jJr 2 2 2 - 6 s 3 (2?3 - C 1C2)] 
S' = 2 ay sc , - ne'e, + ay 


- ac„ + 2ay^s 3 2 3 + -^-8ayj 4 c 2 


a 

t 

a 
b' 

i 

f 

T = 


rco • r x F 


as 2 C2 + 2Sj 3 ?3 + -j-ys 4 c 2 ' 


■ ~ JL ' Qsc i - aj 

= yr ■ 2c„ + aj asc 1 + bs 2 d 2 - Yr 3 (22 3 -CiC2)j 


y£ • 2 ay JC j + ay 


- ac 0 + Ibctjshi + yYDtyj 4 C2 


y L ' ^ c 2 + ay 


as 3 c 3 + -jrbs A c 2 - 27 s 5 (cj- 4S 5 ) 


STEP 4 Numerically Integrate Over Aj To Obtain Elements At s + Aj 


(optional ) 
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STEPS 


s = s + As 
Go back to step 2. 

6.2 Numerical Applications 

The equations of the BG14 5 element method given above in Section 6.1 were programmed as nearly as 
possible in the same format as the older BG14 e method (Bond and Gottlieb, 1989). The two methods 
were then compared to reference cases. The RK45 numerical method (Fehlberg, 1969) was used as the 
numerical integration method in both examples. 

6.2.1 Example j_ 

The first example is that of a highly eccentric ( e ~ 0.95) orbit about the Earth. The orbit is subject to 
the J 2 (Earth oblateness) perturbing potential, which is conservative, plus lunar perturbations. This 
orbit was computed by both BG14 5 and BG14 e methods. This example was also computed by Stiefel 
and Scheifele (1971) with extremely high precision and will be used as the reference. Table I shows 
the components of the position vector in Cartesian coordinates as computed by each method after 50 
revolutions of the satellite. It is seen that both methods compare very closely with the reference but the 
new BG14 5 method being slightly closer to the reference. 

The problem description for the first example is: 

Coordinate system: X and Y fixed in Earth equatorial plane; Z perpendicular to Earth equatorial plane. 


Initial conditions: 


Initial State Vector 

Position 

0.0 


-3400.0 

km 

Velocity 

10.691338 

0.0 

0.0 

km/sec 


The time of comparison is at 288.12768941 days, after approximately 50 revolutions. 


TABLE I - Comparison of BG14 5 and BG14 e Methods 
Final Value Of Position Vector 

Method 

X (km) 

Y (km) 

Z (km) 

Steps/Rev 

(Avg) 

REFERENCE 

Stiefel and Sheifele (1971) 

-24219.0503 

227962.1064 

129753.4424 

500 

BG14 (RK45 Fixed Step) 
8 Method 

-24218.8175 

227961.9146 

129753.3431 

62 

BG14 (RK45 Fixed Step) 
e Method 

-24218.8069 

227961.9186 

129753.3344 

62 


The Earth oblateness and lunar models used are somewhat idealized and are taken from Stiefel and 
Scheifele (1971). These models are specified as follows: 


The Earth oblateness perturbations were compared from the potential model 


3 a t 

V = i(GE)J 2 

2 r J 


Z 2 

2 


3 
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where 


GE = 398601 fan 3 / sec 2 {gravitational constant of Earth) 
a € = 6371,22 km {equatorial radius of Earth) 

J 2 = 1.08265 x 10" 3 ( second harmonic of geopotential) 
The lunar perturbation was computed from 


P = -GM 


r ' - a a 

/r - a/ 3 p 3 


and the lunar ephemeris is approximated by 

a x = p sin O/ 

V3 

a y = — “ p cos Q/ 

= - y p COS Q/ 

p = 384400 km {the Earth-Moon distance) 

Q = 2.665315780887 x 10^ rad /sec {Moon orbital rate) 
GM = 4902.66 km 2 ! sec 2 (gravfrahana/ ca/isfanf a/ Afaan) 


6.2.2 Example 2 

The second example (Adamo, 1989) is that of a near circular geocentric satellite orbit numerically 
integrated by the BG 14 5 method from an initial altitude of 300 km down to entry interface altitude of 
123.278 km (66.565 nautical miles). The perturbations considered were the Jacchia 1970 atmospheric 
model and the GEM- 10 (Lerch, 1979) geopotential restricted to second order and degree. The time of 
flight was about 29.736111 days and the ballistic number was 78.606675 kglm 2 . This case failed at an 
altitude of approximately 135 km (72.894 nautical miles) with the older BG14 e method. 

Coordinate System: True Equator and Greenwich Meridian Of Epoch 


Initial conditions: 


Initial State Vector at UT1 = 0 on 3 September 1991. 

Position 

6677832.962 

-62810.44513 

-27301.63472 

m 

Velocity 

78.98607579 

6821.102837 

3626.863958 

m/sec 


TABLE II - Comparison of BG14 5 and BG14 e Methods 
Final Value Of Position Vector 

Method 

X (m) 

Y(m) 

Z(m) 

Steps/Rev 

(Avg) 

BG14 (RK45 Variable Step) 
5 Method 

2664837.2 

-5838760.8 

1033865.4 

29 

BG14 (RK45 Variable Step) 
e Method 

FAILED 

FAILED 

FAILED 



Additional stress cases (not shown) have been computed in which the solution was propagated down to 
the surface of the Earth (assuming no change in atmospheric density below 90 km). 














7.0 Final Comments 


Recent numerical studies on atmospheric entry from near circular orbits and on low thrust in near circu- 
lar orbits exhibit numerical instability when solved by the method of Bond and Gottlieb (1989) for long 
time intervals. These two cases are similar since both have persistent, tangential, non-conservative per- 
turbations. It was found that this instability was due to secular terms which appear on the right hand 
sides of the differential equations of some of the elements. In this paper this instability is removed by 
the introduction of another vector integral of the motion and another scalar integral which remove the 
secular terms. The introduction of these new integrals require a new derivation of the differential equa- 
tions for most of the elements. For this rederivation the Lagrange method of variation of parameters is 
used making the development more concise. 
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Appendix A - The Variation Of Parameters Method Of Lagrange 


Assume that we have a mechanical system given by 

x=L(x,l) (Al) 

where 

x T = (*i, ' • ' jc n ) 

C = (/V A) 

and t is the independent variable. 


We also assume that the solution of the system of equations (Al) is possible and can be 
expressed 

x = x_{c_ y t) (A2) 

where the integration constants, or parameters, are given by 

cj = (c 1. ■ • ■ ,C„) (A3) 

Now consider another system similar to the system (Al), 

£ -r<£.o + *<2.o ( A4 ) 

where the new term is called a perturbation and is given by 

K t (x, 0 = (gj. • ■ • >g») ( A5 ) 

The objective is to make the solution, equation (A2) of the system (Al), valid for the perturbed 
system (A4) by allowing the parameter £ to be a function of the independent variable. In other 
words the solution (A2) still applies but with the constant (c) replaced by the function ( c(r))* 
So we have 


x = x(c(0» 0 

Now take the total derivative of equation (A 6) 

3x . dx 

X = — c + — 

_ 3c ” dt 

Also take the total derivative of (A2) and use (A 1) to obtain 

f-- ;=/(*.<> 

Note we have used the fact that for unperturbed case the total and partial derivatives of x are the 
same. Using equation (A8) we can eliminate the partial derivative — from equation (A7) 
obtaining. 


(A6) 

(Al) 

(A8) 


£ = | “£ +£(£.0 

Now compare equation (A9) with equation (A4) to obtain 

dx . 


(A9) 


£ 3c 


After the obvious cancellation 


£ + /(£> 0 =£<£» 0 + g(£. t) 


dx . 


(A 10) 


nr 

where the matrix — is obtained from the solution, equation (A2). The matrix must be inverti- 
dc 

ble. That is 
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DET\ 


dx 

Be 


*0 


The system of differential equations for the parameter c is therefore 

ldx W 
c = 


dc 

< —j 


(All) 


282 



Appendix B - The Stumpff Functions 


These functions are related to the trigonometric and hyperbolic functions. The general equation 
for the n th Stumpff function is, 

=,!<-»* (2*T^'" =WA -~ (E1) 

When these series are compared to the series of the trigonometric and hyperbolic functions, the 
following relations exist: 


c 0 {x 2 ) = COS X 

, or 

(-* 2 ) = 

cosh x 

Cj ( x 2) = Sil£ 

, or 

c i(-x 2 ) = 

sinh x 
x 

2v 1 “ COS JC 



cosh x - 1 

)= 2 
jr 

, or 

II 

p' 

r 

x 2 

2 x - sin x 

■ / ” 3 

JC 

r i 


/ 2\ 

sinh x - jc 

, or 

c 3 (~* ) = 

x 3 

r 


cos x - 


c 4 (x 2 ) = 


«-! 


cosh x - 


— , or c 4 (-x 2 ) = 


, * 2 

1 + T 


(B2) 


etc. 


The following identities may also be easily verified: 
c 0 (z) 2 + zcj(z) 2 = 1 

c 0 (z ) 2 - zc x {zf = c 0 (Az) 

c a (z) 2 = 1 - 2zc 2 (4z) (B3) 

c,(z) = 2c 2 (4z) 
c 0 ( z )c i(z) = c i(4z ) 

c 2 (z) = Cj(z) 2 - c 2 (z)c 0 (z) 

The more general identities 

c n+ 2(z) = — Cn+l(z) + — 
n z 

and 

c n (z) + zc n+2 (z) = -y (B5) 

are also valid. 

The derivatives of these functions may be expressed as 

2z ~~T ~ ) = c «-i( z ) “ nc H (z) , n>0 (B6) 

az 

and 

^- = }[«c B+2 (z)-c n+1 (z)] (B7) 

A convenient integration formula is 

Js*c*(ps 2 )ds = s t+I c* +1 (p5^ (B8) 


c B (z) + — c„_i(z) 
n 


, n > 0 


(B4) 
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